function phat = lsidentify(y,u,yl,ul,Q)

% Die Relation y(k) = y(k-1)*a(1)+y(k-2)*a(2)+...+y(k-yl)*a(yl) +
% u(k)*b(0)+k(k-1)*b(1)+...+u(k-ul+1)*b(ul-1) + moeglicherweise
% eine Stoerung v(k) bestimmt fuer jedes k eine Gleichung, die 
% linear in [a,b] ist. Wir wollen das entsprechende 
% (ueberbestimmte) lineare Gleichungssystem loesen um [a,b] zu 
% finden.

firstentry=max(yl+1,ul); % Welches k gibt uns die erste lineare Gleichung
N=length(y);

% S so zusammengestueckelt, dass gilt S*[a.';b.']=ys

ys=y(firstentry:N); 
S = zeros(N-firstentry+1,yl+ul);

for k=1:yl
    S(:,k)=y(firstentry-k:N-k);
end
for k=0:(ul-1)
    S(:,k+1+yl)=u(firstentry-k:N-k);
end

phat = linsolve((S.' * Q(firstentry:end,firstentry:end) * S), (S.'* Q(firstentry:end,firstentry:end) *ys));

end